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ABSTRACT 



On 



We have made a detailed comparison of the results of N-body simulations with the 
analytical description of the merging histories of dark matter halos presented in Lacey & 
^> Cole (1993), which is based on an extension of the Press-Schechter method (Bond et al. 

1991; Bower 1991). We find the analytical predictions for the halo mass function, merger 
rates and formation times to be remarkably accurate. The N-body simulations used 128^ 
particles and were of self-similar clustering, with SI = 1 and initial power spectra P{k) cx 
k" , with spectral indices n = —2,-1,0. The analytical model is however expected to 
apply for abitrary SI and more general power spectra. Dark matter halos were identified 
in the simulations using two different methods and at a range of overdensities. For halos 
selected at mean overdensities ~ 100 — 200, the analytical mass function was found to 
provide a good fit to the simulations with a collapse threshold close to that predicted 
by the spherical collapse model, with a typical error of <^ 30% over a range of 10^ in 
^-^ mass, which is the full dynamic range of our N-body simulations. This was insensitive 

to the type of filtering used. Over a range of 10^ — 10^ in mass, there was also good 
agreement with the analytical predictions for merger rates, including their dependence 
on the masses of the two halos involved and the time interval being considered, and for 
formation times, including the dependence on halo mass and formation epoch. 
The analytical Press-Schechter mass function and its extension to halo lifetimes and 
^ merger rates thus provide a very useful description of the growth of dark matter halos 

through hierarchical clustering, and should provide a valuable tool in studies of the 
C/3 formation and evolution of galaxies and galaxy clusters. 
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1 INTRODUCTION 

In the standard cosmological picture, the mass density of 
the universe is dominated by coUisionless dark matter, and 
structure in this component forms by hierarchical gravita- 
tional clustering starting from low amplitude seed fluctua- 
tions, with smaller objects collapsing first, and then merg- 
ing together to form larger and larger objects. The halos of 
dark matter formed in this way, objects which are in approx- 
imate dynamical equilibrium, form the gravitational poten- 
tial wells in which gas collects and forms stars to produce 
visible galaxies. Subsequent merging of the dark halos leads 
to formation of groups and clusters of galaxies bound to- 
gether by a common dark halo, and is accompanied by some 
merging of the visible galaxies with each other. It is obvi- 



ously of great interest to understand this process of structure 
formation via merging in more detail. One approach, begun 
by Aarseth, Gott & Turner (1979) and Efstathiou & East- 
wood (1981), is to calculate the non-linear evolution of the 
dark matter numerically, using large N-body simulations. A 
second approach, complementary to the first, is to develop 
approximate analytical descriptions which relate non-linear 
properties such as the mass distribution and merging prob- 
abilities of collapsed objects to the initial spectrum of linear 
density fluctuations from which they grew. These analytical 
descriptions must then be tested against the results of the 
numerical simulations. If they work, they provide insight 
into the numerical results, and provide the basis for sim- 
plified calculations and modelling which can cover a much 
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wider range of parameter space than is feasible with the 
numerical simulations alone. 

The analytical approach was pioneered by Press & 
Schechter (1974) (hereafter PS), who derived, rather heuris- 
tically, an expression for the mass spectrum of collapsed, 
virialized objects resulting, via dissipationless gravitational 
clustering of initially cold matter, from initial density fluctu- 
ations obeying Gaussian statistics. The basis of the method 
is to derive a threshold value of the linear overdensity for 
collapse of spherical perturbations, and then calculate the 
fraction of mass in the linear density field that is above 
this threshold when smoothed on various scales. The PS 
mass function formula has since been widely applied to a 
variety of problems, including gravitational lensing by dark 
halos (Narayan & White 1988), abundance of clusters and 
their influence on the cosmic microwave background via the 
Sunyaev-Zel'dovich effect (Cole & Kaiser 1988), and galaxy 
formation (Cole & Kaiser 1989; White & Frenk 1991). An 
alternative derivation of the PS mass function was presented 
by Bond et al. (1991) (hereafter BCEK), who, by consider- 
ing the random walk of linear overdensity (at a fixed loca- 
tion) as a function of smoothing scale, obtained a rigorous 
solution to the problem of above-threshold regions lying in- 
side other above-threshold regions (the so-called "cloud-in- 
cloud" problem). BCEK also showed how the derived mass 
function depended on the filter used to define the spatial 
smoothing, and that the standard PS formula in fact only re- 
sults in the case of "sharp A;-space" filtering, which is also the 
only case for which exact analytical results can be obtained. 
(Related approximate results were also obtained by Peacock 
& Heavens (1990).) In addition, BCEK showed how to go 
beyond a calculation of the mass function at a single time, 
to derive the conditional mass function relating the halo 
masses at two different times. Indepedently, Bower (1991) 
extended the original method used by Press and Schechter 
and derived an identical expression for the conditional mass 
function. This extension to the PS theory was then taken 
up by Lacey & Cole (1993) (hereafter Paper 1), who used 
the conditional mass function to derive a range of results 
on the merging of dark halos. These included the instan- 
taneous merger rate as a function of the masses of both 
halos involved, and the distribution of formation and sur- 
vival times of halos of a given mass identified at a given 
epoch, the formation time being defined as the earlier epoch 
when the halo mass was only half of that at the identifi- 
cation epoch, and the survival time as when the halo mass 
has grown to twice that at the identification epoch. The ex- 
pressions for these depend on the initial linear fluctuations 
through their power spectrum, and on the background cos- 
mology (density parameter Q and cosmological constant A). 
Preliminary applications of these results were made to con- 
straining the value of Q from the merging of galaxy clusters, 
and to estimating the rate of accretion of satellites by the 
disks of spiral galaxies. 

Kauffmann & White (1993) have extended the utility 
of the analytical results of BCEK and Bower (1991) by pre- 
senting a Monte Carlo method of generating merger trees, 
decribing the formation history of halos, that are consistent 
with the analytical conditional mass function. This tech- 
nique has subsequently been utilized in studies of galaxy 
formation by Kauffmann, White & Guiderdoni (1993) and 
Kauffmann, Guiderdoni & White (1994). An alternative 



Monte Carlo implementation has also been used in studies 
of galaxy formation in Cole et al. (1994). 

An alternative analytical approach is to assume that 
objects form from peaks in the initial density field. This has 
been extensively used to study clustering of galaxies and 
clusters (Peacock & Heavens 1985; Bardeen et al. 1986), but 
has been less used in calculating mass functions because of 
the problem of dealing properly with peaks which lie inside 
other peaks, and of identifying what mass object forms from 
a given size peak {e.g. Bond (1988)). Bond & Myers (1993a) 
have recently developed a new method which combines as- 
pects of the Press-Schechter and peaks methods, but this 
requires one to generate and analyze Monte Carlo realiza- 
tions of the linear density field, and so is much more compli- 
cated to apply, though still simpler than doing an N-body 
simulation. 

Our aim in this paper is to test the analytical results 
for merging of dark halos derived in Paper 1 against a set of 
large N-body simulations. Specifically, we test the formulae 
for merger probabilities as a function of the masses of the 
halos involved and of the difference in cosmic epochs, and for 
the distribution of formation epochs of halos. We also revisit 
the question of how well the PS mass function itself works 
compared to simulations. The latter question has been inves- 
tigated previously, notably by Efstathiou et al. (1988) (here- 
after EFWD) for a set of self-similar models with power-law 
initial power spectra and Q = 1, and by Efstathiou & Rees 
(1988), White et al. (1993) and Bond & Myers (1993b) for 
Cold Dark Matter (CDM) models, who all found reason- 
ably good agreement. BCEK tested the conditional mass 
function formulae against EFWD's simulations, and found 
encouraging results, but the size of their simulations (32'^ 
particles) did not allow very detailed comparisons. In this 
paper, we address these questions using 128'^ Si 2 x 10^ par- 
ticle N-body simulations of self-similar models, with Q = 1 
and initial power spectra P{k) oc A;", n = —2,-1,0. With 
this large number of particles, it is possible to make fairly 
detailed tests of the merger formulae. There has been some 
discussion of what is the best filter to use with the standard 
PS mass function formula, and whether improved results 
can be obtained by choosing a threshold linear overdensity 
different from that of the simple spherical collapse model 
(Efstathiou & Rees 1988; Carlberg & Couchman 1989), and 
we investigate this also. Kauffmann & White (1993) made 
a qualitative comparison of some of the properties of the 
merger histories from their Monte Carlo method with those 
from a CDM simulation, and argued that there was reason- 
able agreement. 

An important question that arises when comparing an- 
alytical predictions for mass functions, merger rates etc. of 
dark halos with the corresponding quantities in N-body sim- 
ulations, is how best to identify the halos in the simula- 
tions. Given that the halos found in simulations are nei- 
ther completely isolated nor exactly in virial equilibrium, 
there seems no unique way to do this. Much work has been 
based on the percolation or "friends-of-friends" algorithm 
(Davis et al. 1985) (DEFW), in which particles are linked 
together with other particles into a group if the distance 
to the nearest group member is less than a certain fraction 
(usually taken to be 0.2) of the mean interparticle separa- 
tion. However, other group-finding schemes have also been 
used (e.g. Warren et al. (1992), Bond & Myers (1993b)). It 
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is important to know how much the comparison with analyt- 
ical results depends on the group-finding scheme employed, 
so we will investigate the effect of varying the linking length 
in the friends-of-friends scheme, and also use an alternative 
method based on finding spheres of particles of a certain 
overdensity. 

The plan of the paper is as follows: In §2 we review 
the analytical results on merging derived in Paper 1. In §3 
we describe the N-body simulations, and in §4 we describe 
our group-finding schemes. In §5 we compare the N-body 
mass functions with the PS formula. The following two sec- 
tions test the predictions for merging against the simula- 
tions: merger probabilities as a function of mass in §6, and 
formation epochs in §7. We present our conclusions in §8. 



2 REVIEW OF ANALYTICAL MERGER 
RESULTS 

2.1 Spatial Filtering and Random Trajectories 

In this section, we review the analytical results on merg- 
ing of dark halos derived in Paper 1, BCEK and Bower 
(1991). Central to the approach is spatial filtering of the 
linear density field. The initial conditions for structure for- 
mation are specified as a Gaussian random density field 
8{x) = /9(x)/p — 1 having some power spectrum P{k), where 
p is the mean density. This field is smoothed by convolving 
it with spherically-symmetric filters WR{r) of various radii 
R, having Fourier representations WR{k): 



WR{k) 



WR{k) 



Wr{\x\) exp{-ik.x) dx. 



(2.1) 



The variance of the field after smoothing with a filter is 
related to the power spectrum by 

f CO 

(t^{K) = {8^)r= I Wl{k)P{k)4:7rk^ dk (2.2) 
Jo 

We list below the filters that we will be using in this 
paper, and their Fourier representations. We also give the 
"natural volume" Vf that is associated with a filter of radius 
Ri, defined to be the integral of FF_R(r)/FF_R(0) over all space. 

Top Hat (TH): 



k < 1/Rs 
k > 1/Rs 



(2.10) 



(2.11) 



The filters are normalized according to the condition 
WR(r)4:Trr^ dr = WR(k = 0) = 1. The natural volume of 
a filter is thus Vf = l/FF_R(r = 0). (In the case of the sharp 
A;-space filter, the volume integrals are a little ill-defined if 
done in real space, since the integral FF_R(r)47rr^ dr ac- 
tually oscillates around 1 as r ^ oo. If desired, this minor 
problem can be cured by multiplying WR{r) by exp( — ar) 
before doing the integral, and taking the limit a ^ after- 
wards.) The natural mass under a filter is then defined as 
Mf =-pVi- 

When the density fluctuations are small (8 « 1), they 
grow according to linear perturbation theory, 5(x, t) oc D{t), 
where the linear growth factor D{t) depends on the back- 
ground cosmological model; for Q = 1, D(t) oc a(t) oc t^^^, 
a{t) being the cosmic expansion factor. (We are assuming 
that only the growing mode of linear perturbation theory 
is present.) The non-linear evolution can be calculated an- 
alytically for spherical perturbations (e.g. Peebles (1980); 
Paper 1); for a uniform overdense spherical fluctuation, the 
collapse time depends on its initial linear overdensity. It is 
convenient to work in terms of the initial density field ex- 
trapolated according to linear theory to some fixed reference 
epoch to, for which we also take a{to) = 1; from now on, this 
is what we will mean by 8{x). In terms of this extrapolated 
8, a spherical perturbation of mean overdensity 8 collapses 
at time t ii 8 = 8c{t), where, for Q = 1, we have the usual 
result 



,{t) = 8co/a{t) = 8co{to/t) 
8co = 3(127r)^^V20 w 1.69 



2/3 



(2.12) 



(2.3) 



(The generalization of this for Q < 1 is derived in Paper 1; 
note that the second line of equation (2.1) of that paper 
contains a typographical error.) 

The analytical formulae we now present for the halo 
mass functions, conditional mass functions and lifetimes are 
expressed in terms of this threshold 8c{t), and a{M), the 
variance of the smoothed linear density field as a function 
of the smoothing mass. 



2.2 Mass Function 



WR{k) = [sin(A;i?T) - (kRT) cos(A;i?T)] (2.4) 

(kRT) 



Vt = (47r/3)i?T 
Gaussian ( G): 



WR{r 



1 



WR{k) = exp 



{2^r/^Rl 
k' 



exp 



2RI 



Vg = (2T:fl^R% 
Sharp k- space (SK): 
1 r . 



2RI 



WR{r) 



2^2 r3 



(2.5) 

(2.6) 

(2.7) 
(2.8) 



(r/i?s)-(r/i?s)cos(r/i?s)] (2.9) 



The fraction of mass in halos with mass M , at time t, per 
interval dM , originally derived by Press & Schechter (1974), 
is (Paper 1 equation (2.10)) 



«c(i) 



(27r)i/2a3(M) 



d(T^(M) 



dM 



exp 



^c(t)^ 

2a{Mf 



(2.13) 



a(M) is assumed to decline monotonically with increasing 
M . (To make the correspondence with the equations in Pa- 
per 1, note that we there used the notation S = a^{M), 
Lo = 8c{t).) Thus, the comoving number density of halos of 
mass M present at time t, per dM , is 

dn 



dM 



{M,t) 
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'2^ P S^t) 



din a 



dlnM 



exp 



2(72 (M) 



(2.14) 



where p is now the mean density at the reference epoch to. 
By defining = 8c{t)/a{M), (2.13) can be rewritten as 



d\n 



/ 2 \ 

-= - j/exp(-j/V2) 
J' Vtt/ 



(2.15) 



which is independent of the form of the fluctuation spec- 
trum. 



2.3 Conditional Mass Function and Merger 
Probability 

The conditional probability for a mass element to be part of 
a halo of mass Mi at time ti , given that it is part of a larger 
halo of mass M2 > Mi at a later time t2 > ti, is found by 
considering two different thresholds, 8c{ti) and ^0(^2). The 
result, derived somewhat differently by BCEK and Bower 
(1991), for the conditional mass fraction per interval dMi is 
(Paper 1 equation (2.15)) 



df 

dMi 



{Ml,tl\M2,t2) 

{8ci - 8,2} 



)3/2 



dal 



dMi 



exp 



{8ci 



(2.16) 



where ui = a{Mi), U2 = cr{M2), 8ci = 8c{ti) and 8c2 = 
^0(^2). The only assumption made here about the function 
8c{t) is that it monotonically decreases with increasing t. 
The reverse conditional probability, for M2 given Mi , is (Pa- 
per 1 equation (2.16)) 



df 

dM2 



{M2,t2\Ml,tl 

1 8,2(8c 



(27r)i/2 



3/2 



du\ 




dM2 


exp 



{8c2<jj 



cl'72) 



(2.17) 



This is obviously the same as the probability for a halo of 
mass Ml at ti to be incorporated into a halo of mass M2 > 
Ml at t2 > ti. Thus, if we set M2 = Mi + AM and t2 = 
ti -\- At in the above formula, we get the probability for a 
halo to gain mass AM by merging in time At. Taking the 
limit At then gives the instantaneous merger rate as a 
function of Mi and AM (equation (2.18) of Paper 1). 



2.4 Formation Times 

Suppose one identifies a halo of mass M2 at time t2 . At an 
earlier time, one can identify the progenitors of this halo. We 
define the formation time of the halo identified at epoch t2 as 
the earliest time ti < t2 at which it has a progenitor of mass 
Ml at least half of M2 . We find the cumulative probability 
distribution for ti (equation (2.26) of Paper 1) 

P{tf < tl\M2,t2) = 

l-Mn 



M2 df 



' M2/2 



Ml dMi 



{Mi,ti\M2,t2)dMi 



(2.18) 



The differential probability distribution for ti is then given 
by 

J^(t.|M2,t2) = 



M2/2 



M^ 

Ml 



d r df 



dtf IdMi 



{Ml,ti\M2,t2) 



dMi 



(2.19) 



We noted in Paper 1 that the expression corresponding 
to equation (2.19) actually leads to a slight mathematical 
inconsistency in some cases: for power-law power spectra 
P{k) oc with ra > 0, the probability density for ti goes 
slightly negative for small t2 — tf. In Paper 1, we also de- 
rived formation time distributions based on a Monte Carlo 
method, which do not have the problem of negative proba- 
bility density. These distributions have similar shapes to the 
analytical ones, but with the mean shifted. We will see in §7 
that the analytical distribution gives a remarkably good fit 
to the N-body results. 



2.5 Self-Similar Models 

The analytical results presented above make no special as- 
sumptions about the functions a{M) and 8c{t), except that 
they are monotonic. However, in testing these results against 
simulations, we will focus on self-similar models, in which 
the density parameter Q = 1, so that there is no char- 
acteristic time in the expansion of the universe, and in 
which the initial density fluctuations have a scale-free spec- 
trum, P{k) oc k" . In this case, the evolution of struc- 
ture should be self-similar in time. This has some advan- 
tages, to be discussed in §3. From equation (2.2), we obtain 
a(M) oc M~^"'^^^^^ , in general. For a to decline with in- 
creasing M, we require n > —3, which is just a condition 
for structure to grow hierarchically, with small objects col- 
lapsing first and then merging to form larger objects. For 
the top hat filter, the integral in equation (2.2) only con- 
verges for ra < 1. We will be considering N-body models 
with ra = -2, -1, 0. 

For Q = 1, density fluctuations grow as D{t) oc a{t) oc 
t^^^ in linear theory, so that t he r.m.s. fl uctuation on co- 
moving scale k~ is roug hly y/4:Trk^P{k)a{t). Thus we can 
define a characteristic non-linear wavenumber A;^(t) by 



A7rkl{t)P{k4t))a{tf 



1 



(2.20) 



At time t, fluctuations are of order unity and are starting 
to collapse on a comoving lengthscale ~ k^(t)~^ . In a self- 
similar model, this should be the only characteristic length- 
scale for structure. 

In our analytical expressions for mass functions, merger 
probabilities etc, it is convenient to define a filter-dependent 
characteristic mass scale M^(t) by 



aiM^t)) = 84t) 



(2.21) 



where 8c{t) = 8co / a{t) for Q = 1. The mass is related to 
the filter-independent quantity A;^ by 



M^{t) 



7f 



, 6/(n + 3) _ 

Cf(ra) \ p 



klit) 



(2.22) 



where p is the mean density at the reference epoch when 
a = 1. 7f relates the volume of a filter to its radius through 
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Vf = ffRf (c.f. §2.1), and Cf(ra), which enters through the 
relation (2.2) between <t{R) and P{k), is defined by 



Cf (ra) 



For the 

for TH 
and (2. 

a6/(n + 3) 

as 



(2.23) 



filters we are using, c/(-2) = 1.373,0.941,1.000, 
= 1.500,0.707,0.707 and c/(0) = 2.170,0.666,0.577 
, G, SK filters respectively. From equations (2.20) 
22), the characteristic mass grows as M^(t) oc 
The mass function (2.13) can then be rewritten 



df 



dlnM 



1/2 



ra + 3 



X exp 



M_ 

6 J 

X (n + 3)/3 



(n + 3)/6 



(2.24) 



in which form the mass and time only appear in the com- 
bination M I Mi^{t). Similarly, the merger probability (2.17) 
can be rewritten as a distribution for M2 1 Mi depending on 
Ml I Mi^(ai) and a2/ai, and the formation time distribution 
(2.19) can be rewritten as a distribution for ai I a2 depending 
on M2/M*(a2). 



2.6 Choice of filtering and collapse threshold 

The BCEK derivation of the PS mass function and of the 
conditional mass function is based on sharp A;-space fil- 
tering. On the hand, the original, more heuristic, deriva- 
tion by PS themselves assumed top hat filtering, as did 
Bower's derivation of the conditional mass function. Most 
applications of the PS formula have followed the latter ap- 
proach and used the a(M) relation for top hat filtering (with 
M = (4:w/3)'pR'^), but some have instead used a(M) for 
Gaussian filtering (with M = (2Tr)^ ^^'pR%) (e.g. Efstathiou 
& Rees (1988)). Even for a given choice of filter, one can 
obtain different a{M) relations simply by choosing a differ- 
ent mass-radius relation from the "natural" one discussed in 
§2.1; after all, it is not obvious how the mass of a collapsed 
object is related to the profile of the filter used to identify 
it. BCEK suggested calculating the filter mass for a general 
filter from M = (47r/3)p-RT, where the "eqivalent top hat" 
radius Rt is defined through the relation cr{R) = aTH{RT), 
on the grounds that the collapse threshold 8c{t) is also cal- 
culated for a top hat spherical perturbation. For power-law 
power spectra, this requires making ji in equation (2.22) a 
function of spectral index as well as filter type, while for a 
general P{k) it must be a function of R. This procedure is 
equivalent to using the a{M) relation for top hat filtering 
in formulae like (2.13) and (2.17), even if these formulae are 
derived for sharp A;-space filtering. 

In this paper, we adopt an empirical approach to the 
choice of filter and M{R) relation. We will compare the re- 
sults of N-body simulations to the formulae using top hat, 
Gaussian and sharp A;-space filtering for a{R). We assume 
a mass-radius relation M = fi'pR'^ , with fi a constant de- 
pending on the filter type but independent of the power- 
spectrum, but will consider values of ji different from the 
"natural" ones given in §2.1. 

A related issue concerns the choice of collapse threshold, 
which for Q = 1 boils down to a choice for 8co in equation 
(2.12). While the spherical collapse model gives an unam- 



biguous answer, one can take the view that, since real col- 
lapses have non-top hat and non-spherical density profiles, 
8co should be regarded as a phenomenological parameter, 
chosen to give the best fit to N-body results (e.g. Bond 
& Myers (1993b)). Since the PS mass function and other 
formulae depend on ji and 8co only through (equation 
2.22), there is a degeneracy between these two parameters 
for any given spectral index ra, but this degeneracy is lifted 
as soon as one considers results for different ra. The choice of 
8co and ji will be considered in relation to the simulations 
in §5. 



3 N-BODY SIMULATIONS 

The simulations were performed using the high resolution 
particle-particle-particle-mesh (P^M) code of Efstathiou et 



al. (1985) (EDFW) with 128'' 



2 X 10" particles. The 



long-range force was computed on a 256^ mesh, while the 
softening parameter for the short-range force was chosen to 
he = 0.2(i/256), where L is the size of the (periodic) 
computational box. This corresponds to the interparticle 
force falling to half of its point-mass value at a separation 
r K, 0.4); K, i/3200. Initial positions and velocities were gen- 
erated by displacing particles from a uniform 128'^ grid ac- 
cording to the Zel'dovich approximation, assuming the linear 
power spectrum and Gaussian statistics. We ran one simu- 
lation for each of ra = —2,-1,0. For ra = — 1 and ra = 0, 
the initial amplitude of the power spectrum was chosen to 
equal the white noise level at the Nyquist frequency of the 
particle grid; for ra = —2, this choice was found to lead to 
large departures from self-similar behaviour in the derived 
mass functions and related quantities, so this simulation was 
instead started when the amplitude was smaller by a factor 
0.4. We adopted the convention of normalizing the expansion 
factor a to 1 when the variance of the linear theory power 
spectrum in a top-hat sphere of radius i/32 was 1. With this 
choice, the initial expansion factors were a, = 0.2,0.15,0.06 
respectively for ra = —2, —1, 0. The initial r.m.s. ID displace- 
ments of the particles were approximately 0.9, 0.5, 0.25 in 
units of the particle grid spacing. The time integration was 
performed using the variable p = a" , with a = 2/(ra -|- 3), 
and a constant stepsize l^p/pt = ja(256i]/L)^^^ x 0.023a 
(EDFW). 

For each simulation we output the positions and ve- 
locities of all the particles at many epochs. The spacing of 
these outputs was chosen so that the characteristic mass, 
M^, increased by a factor \/2 between each successive out- 
put, which corresponds to an increase of a factor 2'-""'"'^-'^"'^ 
in the expansion factor a. The final expansion factors for 
the simulations were determined basically by the computer 
resources available; for the later stages, the CPU time was 
dominated by the short-range force calculation by a large 
factor. The simulations were stopped at expansion factors 
a = 1.26, 2.00, 1.68 respectively for ra = —2, —1, 0. This gave 
us between 20 and 24 useful outputs from each simulation. 

Self-similar models have two advantages from the point 
of view of testing our analytical predictions against simula- 
tions, (i) We can check whether the simulations obey the self- 
similar scaling which physically they should. In particular, 
this allows us to test whether the simulations were started 
at a small enough expansion factor. We can also delineate 
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the range of masses of virialized objects for which the sim- 
ulations are reliable, bearing in mind the effects of particle 
discreteness and force resolution on small scales, and missing 
long-wavelength modes on scales larger than the box. (ii) By 
using the self-similar scaling, we can straightforwardly com- 
bine the results from different output times to reduce the 
Poisson fluctuations on the N-body results due to the finite 
number of halos in the box. These factors, and the simple 
way to test for dependence on the form of the initial spec- 
trum, were what motivated us to look at self-similar models 
rather than more physically-inspired models such as Cold 
Dark Matter. 

4 GROUP-FINDING IN THE SIMULATIONS 

4.1 Overview 

We wish to obtain the properties of dark matter halos from 
the simulations to compare with our analytical results. Sim- 
ple theoretical calculations idealize dark halos as isolated 
spherical objects in dynamical equilibrium, but the objects 
found in simulations with Q = f are neither isolated nor in 
complete dynamical equilibrium, because halos continue to 
grow by accreting or merging with other halos on a timescale 
comparable to the expansion time, which is also comparable 
to their internal dynamical times. (In a universe with Q ^ f , 
on the other hand, one expects the conditions of isolation 
and dynamical equilibrium to be much better satisfied). Nor 
are real halos spherical. The issue of how to identify the 
groups of particles in the simulations that one calls dark ha- 
los is therefore not completely straightforward, and a variety 
of schemes have been used by different authors (e.g. DEFW, 
Barnes & Efstathiou (f987). Warren et al. (f992). Bond & 
Myers (f993b)). In order to give some idea of how sensitive 
the comparisons are to the group-finding method employed, 
we will present results for two different schemes, namely 
percolation, and a spherical overdensity method. Both are 
based on particle positions, making no use of the velocity 
information in the simulations. Both methods effectively in- 
volve choosing a density threshold to define groups (a local 
density in one case and a mean density inside a sphere in 
the other), but do not build in any preferred length or mass 
scales. The latter is important if one wishes to study proper- 
ties like the distribution of halo masses. A fuller discussion 
and comparison of the internal properties of the halos found 
by using these different schemes and parameters, for the dif- 
ferent initial power spectra, will be given by Cole & Lacey 
(f994). 

4.2 Friends-of-Friends (FOF) Groups 

The percolation method is the standard friends-of-friends al- 
gorithm (hereafter FOF) of DEFW, which has been widely 
used. Groups are defined by linking together all pairs of 
particles with separation less than b'n~^^^ , W being the mean 
particle density. This defines groups bounded approximately 
by a surface of constant local density, p/Ji Si 3/(2x6'^). The 
value usually used for the dimensionless linking length is 
b = 0.2 (e.g. Frenk et al. (f988)), which corresponds to a 
density threshold p/p ~ 60. For a spherical halo with a den- 
sity profile p(r) oc , this local density threshold corre- 



sponds to a mean overdensity (p) /"p ~ f80, which is close to 
the value f87r^ Si f78 for a just virialized object predicted for 
a top-hat spherical collapse (e.g. Peebles (f980). Paper 1). 
It has been argued that the choice b = 0.2 approximately 
delineates between objects which are virialized and objects 
which are still collapsing in their outer parts, but in our own 
investigations (Cole & Lacey f994), we find that the close- 
ness to global virial equilibrium of N-body halos depends 
rather weakly on the value of b, so that this criterion does 
not strongly select a value for b. We will make comparisons 
with FOF groups for b = 0.15, 0.2, 0.3, corresponding to lo- 
cal overdensities of f40, 60 and f 8 respectively. We will use 
the abbreviation FOF(6) for groups identified using FOF 
with linking parameter b. 



4.3 Spherical Overdensity (SO) Groups 

The second method we apply, which we call spherical over- 
density (SO), is based on finding spherical regions in a sim- 
ulation having a certain mean overdensity, which we denote 
by K = (p) /"p- We first calculate a local density for each 
particle by finding the distance rn to its N'th nearest neigh- 
bour, and define the density as 3(N -|- 1)/ (Awr'^). Particles 
are sorted by density. The highest density particle is taken as 
the candidate centre for the first sphere. A sphere is grown 
around this centre, with the radius being increased until the 
mean overdensity first falls below the value k. (The sphere 
must contain at least 2 particles). The centre of mass of the 
particles in this sphere is then taken as a new centre, and 
the process of growing a sphere of the specified overdensity 
repeated. This process is iterated until the shift in the centre 
between successive iterations falls below c'n~^^^. The parti- 
cles in the sphere are all labeled as belonging to the same 
group, and removed from the list of particles considered by 
the group-finder. Then one moves on the next highest den- 
sity particle which is not already in a group, and repeats 
the process of finding an overdense sphere, including iter- 
ation of the centre. Finally, after all the groups have been 
found, any groups which lie inside larger groups are merged 
with the larger group, according to the following procedure: 
Each group is considered as a sphere of mean overdensity k, 
based on its actual mass, centred on its actual centre of mass. 
Starting with the largest sphere, any smaller sphere whose 
centre is inside the larger sphere is merged with that sphere 
(but the assumed radius of the larger sphere is not changed). 
This is then repeated for the next largest remaining sphere, 
and so on down in mass. We will use the abbreviation SO(k) 
for SO groups identified at overdensity k. 

We will present comparisons for SO groups with spher- 
ical overdensity k = f80, chosen to agree with the spheri- 
cal collapse model. For the local density estimate, we used 
N = fO, but the results were not especially sensitive to this. 
For the convergence of the group centres, we found e = O.f 
to be adequate. With these parameters, an average of fewer 
than O.f iterations per group were required, and the frac- 
tion of particles involved in merging small groups into large 
ones was less than fO~'^. We also experimented with defining 
the initial list of centres for growing spheres from particles 
ranked either by gravitational potential (deepest potential 
first) or by the magnitude of the acceleration (highest ac- 
celeration first). The potential and acceleration were calcu- 



Merger rates 7 



lated for each particle using a modified version of the P"^ M 
code, with the same grid and softening parameters as for 
the original simulation. The groups found starting from ac- 
celeration or potential centres were almost identical to those 
found starting from density centres, and the mass functions 
agreed to within a few percent. A fuller discussion of the 
method will be given in Cole & Lacey (f994). We have used 
the method based on density centres in this paper since it 
is more straightforward. 

Our spherical overdensity algorithm has many similar- 
ities to the method used by Warren et al. (f992), and to 
the "smooth particle overdensity" method of Bond & Myers 
(f993b). The Warren et al. method grows spheres of speci- 
fied overdensity around centres which are particle positions 
ranked by gravitational acceleration, and merges small ha- 
los which are inside larger ones, but does not iterate the 
positions of the centres. The Bond & Myers method grows 
spheres of a given overdensity around centres which are cho- 
sen initially to be positions of particles ranked by local den- 
sity, and then iterates each sphere until the centre of mass 
converges. However, the volume of a group, used in defining 
its mean density, is defined by calculating a volume for each 
particle based on its local density, and summing these over 
all particles within the sphere to get the total volume, rather 
than as the volume of the sphere itself. 



4.4 Comparison of Methods 

The percolation method has the advantage that it is simple, 
relatively fast, and does not make any assumption about the 
geometry of the groups. However, some (a definite minor- 
ity) of the groups it selects are formed of two or more dense 
lumps separated by low-density bridges, as has been noted 
by previous authors. These groups seem rather unphysical. 
The spherical overdensity method avoids this problem, in- 
stead producing groups concentrated around a single centre, 
but on the other hand, tends to chop off the outer portions 
of ellipsoidal halos. It is also more complicated and more 
time-consuming to run. It is not clear which method should 
be considered "best" , so it is interesting to compare results 
obtained with both. 



5 MASS FUNCTIONS IN THE SIMULATIONS 
5.1 Tests for Self- Similarity 

Figure f shows the N-body mass functions for all output 
times for all three values of the spectral index n. The halos 
were found using the friends-of- friends method with h = 0.2, 
and the mass function for each output time rescaled to a dis- 
tribution in MIMi,(t), with M*(t) oc a(tfl^"+^'> computed 
using a top hat filter with the standard choices n = 47r/3 
and 8co = 1.69. Scaled in this way, the mass functions 
for different times should be identical, as a simple conse- 
quence of self-similarity, and it can be seen that the sim- 
ulations conform to this expectation very well. We have 
excluded from the plots halos having N < iVmin = 20 or 
N > iVmax = 2 X 10* particles. Halos containing only a 
few particles are not represented accurately by the simula- 
tions, and cause noticeable departures from self-similarity 
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Figure 1. Mass functions for different output times, for groups 
identified using FOF(0.2), and masses rescaled to be in units of 
the characteristic mass M^[t). df fdlnM is the fraction of mass 
in halos per logarithmic interval in halo mass. Each panel shows 
the results for a single N-body run (ri = — 2, — 1, 0), with the mass 
functions for different output times plotted as solid lines, and the 
Press-Schechter prediction (equation (2.24)) for top-hat filtering 
and Sco = 1.69 as a dashed line. 
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if included, while the properties of halos having masses ap- 
proaching that of the box are affected by the absence from 
the initial conditions of modes with wavelength exceeding 
the size of the box. For our simulations, the iVmax cutoff in 
fact makes little difference, because is still much smaller 
than the mass of the box when the simulations are halted. 

Rigorously, self-similar clustering solutions only exist 
for spectral indices in the range — f < n < f, since out- 
side this range, the r.m.s. peculiar velocity receives diver- 
gent contributions from either large or small scales (Davis 
& Peebles f977). For n < — f, the problem arises from the 
very long wavelength fluctuations, which contribute negligi- 
bly to the r.m.s. overdensity, but produce large-scale coher- 
ent velocities. This should not affect the collapse of structure 
on small scales, so small scale structure is still expected to 
evolve in a self-similar way. This expectation is borne out 
by the scaling of the mass functions, which works as well for 
the n = —2 models as for ra = — f or 0. 

Since the scaling of the mass functions in the simula- 
tions is consistent with self-similarity, we average the re- 
sults for all output times together to increase their preci- 
sion, weighting the contribution to df/dlnM in each bin in 
M I Mi^ in proportion to the number of halos in that bin for 
each output time. We use these averaged mass functions in 
the comparisons that follow. 



5.2 Choice of Filter and Density Threshold 

We now consider the choice of filter and collapse threshold in 
the Press-Schechter mass function (c,/, §2.6). The PS pre- 
diction is that when the distribution in mass is converted 
to a distribution in p = 8c{t)/a{M), it should have a uni- 
versal form (equation 2.f5), independent of the power spec- 
trum. However, a{M), and thus j/, depend on the filter used. 
For self-similar models, j/ = (M/M^(^))("+^)/^ and differ- 
ent choices of filter lead to different values of M^, according 
to equation (2.22). For a given filter, the value of also 
depends on the parameters 8co and ji (recall that the filter 
mass is related to the radius by Mf = fiJiRy). 

We first test how well the N-body mass functions for dif- 
ferent initial power spectra agree with each other when ex- 
pressed as functions of j/. Figure 2 shows the mass functions 
of FOF(0.2) groups in our self similar models, converted to 
distributions in j/, for three different choices of filter: top 
hat (TH), Gaussian (G) or sharp A;-space (SK). The relative 
placement of the N-body curves for different n depends on 
7f, but is independent of 8co, which just causes a uniform 
shift of all the curves by the same amount. For top hat and 
sharp A;-space filtering we have used the "natural" values of 
7f from §2.f , while for Gaussian filtering, we have used both 
the natural value and one 2.5 times larger. We have assumed 
8co = f .69. It can be seen that for top hat and sharp A;-space 
filtering, the N-body mass functions df / din v for different n 
agree fairly well using the natural values for 7f (47r/3 and 
67r^ respectively), so we have not considered different values. 
For Gaussian filtering, on the other hand, there is a large 
spread between the N-body curves for different n when 7f is 
taken to be its natural value (27r)'^^^, but good agreement if 
the filter mass is increased by a factor 2.5. Previous authors 
who have used Gaussian filtering in conjuction with the PS 
formula {e.g. Efstathiou & Rees (f988), Carlberg & Couch- 



man (f989)) all appear to have assumed 7f = (27r)'^^^. The 
effect of the latter choice is that the best fitting value for 8ca 
depends on the shape of the power spectrum. These results 
about the best values for 7f for different filters were found to 
apply equally well for the other group identification schemes 
we tested (friends-of- friends with h = 0.15 and b = 0.3, and 
spherical overdensity with k = fSO). 

Having chosen optimal values for ji, a second question 
is: how well does the PS formula actually fit the N-body 
results, and what is the best value to use for 5cO? In Fig- 
ure 2, which assumes 8co = 1.69, the standard PS prediction 
(equation 2.f5) is shown by a dot-dash curve in each panel. 
The PS and N-body curves have very similar shapes, but 
it is apparent that they could be brought into even closer 
agreement by shifting the N-body curves horizontally, which 
corresponds to changing 8co (i' 8co). For each filter, we 
have estimated by eye what value of 8co gives the best fit of 
PS with N-body mass functions. The best fitting PS curve 
is shown by a dotted line in each panel, and the correspond- 
ing value of 8co displayed. (Rather than replot the N-body 
curves with the new value of 8co, we have shifted the PS 
curve inversely as, u oc 8~g .) For Gaussian filtering with the 
non-optimal ~/f = (27r)'^^^, the best fit 8co depends strongly 
on the power spectrum, so we have fitted the n = —1 results, 
in the middle of our range. Using the optimal ji values, it 
can be seen that the best fit 8co values for FOF(0.2) groups 
differ by less than 20% from the canonical 8co = f.69 for 
each of the three filters considered. When these best fit val- 
ues are used, the PS formula fits the N-body mass functions 
extremely well, with an accuracy better than 30% over the 
range 0.3 < i' < 3. However, in all cases the PS formula sys- 
tematically over-estimates the abundance of low mass halos 
(j' < f) compared to the simulations. 

These results on 8co are reasonably consistent with 
those found by previous authors. EFWD found a reason- 
able fit to FOF(0.2) groups in self-similar models for top 
hat filtering with 8co = 1.68, while Efstathiou & Rees 
(f988) and Carlberg & Couchman (f989) found 8co = f.33 
and 8co = f.44 respectively for FOF(0.2) groups in CDM 
models using Gaussian filtering with ~/f = (2Tr)^^^ (note 
that the CDM power spectrum has effective spectral in- 
dex din P{k)/dln k Si —2 on the range of scales consid- 
ered). Bond & Myers (f993b) find 8co = f.58 from their 
CDM simulations, but this is using their own group-finding 
scheme, which produces somewhat more massive groups 
than FOF(0.2), and correspondingly requires a smaller 8co. 



5.3 Mass Functions with Different 
Group-finding 

It is important to investigate how sensitive the results on 
mass functions are to the group-finding method employed. 
We have repeated the previous analysis for friends-of-friends 
with two other values of the linking parameter, b = 0.15 and 
0.3, and for the spherical overdensity method with overden- 
sity K = f 80. Figure 3 shows the results for top hat filtering 
with 7f = 47r/3 and 8co = 1.69. The dot-dash curve in each 
panel is identical, and is the standard PS result, while the 
dotted curve is the shifted PS curve which seems to give the 
best fit, the corresponding best fit value of 8co being given 
in each panel. Comparing the N-body mass functions with 
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Figure 2. Averaged mass functions for the different N-body models, compared to the Press-Schechter mass function using different 
fihers. df /dlni' is the fraction of mass per logarithmic interval in = (M /M^Y'^'^'^^ . Groups were identified using FOF(0.2). Each 
panel shows results for a different choice of filter: top hat (TH), Gaussian (G) and sharp A;-space (SK) respectively. For Gaussian filtering, 
we show results for two different values of 7f , which relates the filter mass to filter radius through Mj^ = 7f"p/?^. Each panel shows the 
N-body mass functions for n = —2 (long dashed line), n = —1 (solid) and n = (short dashed), where these have been averaged over all 
output times. Also shown is the standard Press-Schechter result (equation (2.15)) (long dash dot curve), and the PS curve shifted in v 
to give the best fit to the N-body results (dotted curve). Shifting in v is equivalent to using a different value for Sco from the standard 
one. The values of 7f assumed, and of Sco for the best-fitting PS curves, are shown in the lower left corner of each plot. 
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Figure 3. Effects on N-body mass functions of different group-finding scliemes. Eacli panei sliows tlie mass functions for n = —2, —1, 
(long daslied, solid and short dashed lines respectively), for a different choice of group finding scheme and parameters. The first 3 panels 
show friends-of-friends with b = 0.15,0.2,0.3 respectively, while the last shows the spherical overdensity method with k = 180. In each 
panel, the N-body mass functions have been averaged over output times, and converted to distributions in v assuming top hat filtering 
with 7f = 47r/3 and Sco = 1.69. Also shown in each panel is the standard Press-Schechter result (equation (2.15)) (long dash dot curve), 
and the PS curve shifted in v to give the best fit to the N-body results (dotted curve). The values of Sco corresponding to the best-fitting 
PS curves are shown in the lower left corner of each plot. 
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different group finders, tfie piots sfiow tfie expected resuit 
tfiat fiaio masses are smaiier for FOF(0.15) tfian FOF(0.2), 
and iarger for FOF(0.3). SO(180) fiaios are aiso on average 
smaiier tiian FOF(0.2), but give very simiiar mass functions 
to FOF(0.15). Tlie best fit vaiues of 8co vary witli group- 
finder accordingiy, iarger liaio masses requiring a smaiier 
8co, and vice versa, as siiown in tlie figure. Witli TH fiiter- 
ing, aii tlie group-finders give reasonable fits to PS (witli 
different ^co), though in each case, when the PS formula is 
fitted at the high mass end, it predicts somewhat too many 
low mass halos, the discrepancy being largest for FOF(0.15) 
and SO(180). FOF(0.3) groups result in the best fit to PS 
overall, when 8co is allowed to vary. For the canonical choice 
8co = 1.69, FOF(0.2) and FOF(0.3) groups agree about 
equally well with PS. For Gaussian and sharp A;-space fil- 
tering, with 7f = 2.5 (27r)'^^^ and Gtt^ respectively, the fits of 
the PS formula to the N-body mass functions for the differ- 
ent group finders are about as good as for top hat filtering 
when 8co is allowed to vary. 

5.4 Discussion 

To summarize the results of this section: at least for power- 
law initial power spectra, P{k) oc with —2 < ra < 0, there 
seems little to choose between top hat, Gaussian or sharp 
A;-space filtering in the Press-Schechter mass function, pro- 
vided the factor ji is chosen appropriately. For top hat and 
sharp A;-space filtering, the natural values (47r/3 and Gtt^ re- 
spectively) work well, while for Gaussian filtering, a value 
around 2.5 (27r)'^^^ seems best. The N-body mass functions 
change with changes in the group-finding method or in its 
parameters, but this can to a considerable extent be com- 
pensated for in the PS formula by adjusting 8co- If ^co is 
instead fixed at the value 1.69 from the spherical collapse 
model, then for FOF groups, a value b Si 0.2 — 0.3 seems to 
give the best agreement with PS. We caution however that 
these conclusions may be changed for more general power 
spectra, in particular, there may be more difference between 
the top hat filter, which cuts off only as Wn^k) oc k~^ at 
short wavelengths, and the Gaussian or sharp A;-space filters, 
which both cut off much more sharply. In the remainder of 
this paper, in looking at merger properties, we will concen- 
trate on 6 = 0.2 FOF groups, using top hat filtering with 
7f = 47r/3, since these choices are the most standard. For 
simplicity we will use the canonical value 8co = 1.69 for the 
collapse threshold, since this is in any case close to the best 
fit for the other choices made. The PS mass function then 
fits our N-body results to within a factor 2 or better for 
0.3^1/^3. 



6 MERGER RATES 

The conditional mass function, df /din AM\Mi,Aa,ai, (2.17) 
describes the probability that a halo of mass Mi at the 
epoch when the expansion factor equals ai will accrete mass 
AM = M2 — Ml to become a halo of mass M2 at expansion 
factor a2 = ai -\- Aa. In the limit of Aa this yields the 
instantaneous merger rate at expansion factor ai of halos of 
mass Ml with halos of mass AM. Here we compare the con- 
ditional mass function (2.17) for both small and large time 



intervals At with estimates made directly from our N-body 
simulations. 

Having constructed group catalogues for each output 
time of our simulations using one of the group finding al- 
gorithms discussed in §4 it is an easy matter to construct 
the conditional mass function for any pair of output times. 
For each group of mass Mi identified at the epoch when the 
expansion factor equals ai we determine which halo it has 
become incorporated in at expansion factor a2 by finding the 
halo at this later epoch which contains the largest fraction 
of the particles from the original halo. Defining the mass of 
this new halo as M2 we construct a joint histogram of Mi 
and AM = M2 — Mi . The scale free nature of the initial 
conditions of our simulations implies that the form of these 
histograms when expressed in units of Mi /M^ (Af^ given 
by (2.21) at a = ai) and AM I Mi should depend only on 
Aa/ai = (a2 — ai)/ai. Consequently the histograms from 
different pairs of output times, but with the same value of 
Aa/ai, can be averaged to yield more accurate estimates 
of the conditional mass function. For each pair of output 
times we weight the contribution to df /din AM\Mi,Aa,ai in 
each bin of AM /Mi in proportion to the number of halos in 
that bin. We estimate Poisson error bars for the conditional 
mass function from the number of halos Nhin in each mass 
bin, summed over the pairs of output times. Successive pairs 
of output times are not really independent, so we define an 
effective number per bin as Nes = Nhin/ f , and take the frac- 
tional error to be 1 /V NeS- For the particular case for which 
Aa/ai corresponds to the spacing between successive out- 
puts, i.e. outputs spaced by a factor t/2 in M^(t), we take 
/ = 1, but for all larger values of Aa/ai we adopt / = 2, 
which corresponds to taking outputs separated by a factor 
2 in M^(t) to be independent. This procedure is obviously 
not rigorous, but provides some indication of the magnitude 
of statistical errors. 

We have compared the analytical predictions of (2.17) 
with the conditional mass functions estimated from each 
of our the simulations for a wide range of both Mi and 
Aa/ai. A representative selection of these comparisons are 
shown in Figures 4, 5 and 6, which are for the case of groups 
identified using FOF(0.2), and restricted to halos satisfying 
Ni > 20 and AN = N2 - Ni > 20, Ni and N2 being the 
number of particles in the halo at ai and a2 respectively. In 
these figures, was defined using top hat (TH) smoothing 
with 7f = 47r/3 and the standard 8co = 1.69, which gave 
close to the best fit to the FOF(0.2) mass functions in §5. 
The fits of the analytical to the N-body results in these 
diagrams are much less sensitive to the adopted value of 
8co than was the case for the mass functions, and adopting 
the "best fit" value of 8co = 1.81 only slightly changes the 
analytic distributions. In fact, the values of 8co which give 
the best fit for the conditional mass functions are in general 
somewhat different from those found for the mass functions 
themselves. The deviations of the first few or last few N- 
body points from the analytical curves seen in some of the 
plots seem not to be significant, but depend on the choice of 
Ni and AN. With more conservative cuts, one gets smaller 
deviations, but then the N-body data covers smaller ranges 
in Ml and AM. 

The conditional mass functions, df /din AM\Mi,Aa,ai, 
shown in these figures exhibit a quite complicated depen- 
dence on the accreted mass, AM, the initial mass. Mi, the 
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Figure 4. A comparison of the analytical conditional mass functions (2.17) (dashed curves) with those of FOF(0.2) groups in the n = —2 
simulation (solid lines). Top hat (TH) smoothing and the standard Sco = 1.69 were used to define M^. We have used only halos with 
A^l > 20 and AA^ > 20. An indication of the magnitudes of the statistical errors in these estimates are shown by the Poisson error bars, 
which are calculated as described in the text. The top row of plots show results for an interval Aa/ ai = 0.06, which equals that between 
consecutive outputs of our n = —2 simulation, i.e. between which increases by a factor of \/2- From left to right these three plots 
correspond to increasing Mi /M^ as indicated on each plot. The lower row of plots make the same comparison for a larger time interval, 
Aa/ai = 0.59, which corresponds to the interval over which increases by a factor 16. 
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Figure 5. Like Figure 4, but for the n = —1 model. The values of Aa/ ai indicated on the figures again correspond to the intervals over 
which increases by factors \/2 and 16 respectively. 



time interval, Aa, and on the spectral index, n, of the ini- 
tial conditions. The distributions in AM/Mi are asymmetric 
and vary in both the form and degree of this asymmetry, in 
their width, and in the location of their peak, when either 
Ml / , Aa/ai, or n are changed. The distributions are 
broadest and most asymmetric for low values of Mi/M^. 
Also for fixed Mi /M^ they are broader for lower values of 
n. Increasing Aa/ai shifts the peaks of the distributions to 
higher masses while also altering the asymmetric nature of 
the low Mi/Mi^ distributions. Remarkably all these features 
are quantitatively reproduced by the analytical expression 
(2.17). Overall, we find reasonable agreement between the 
simulations and the analytical distributions over 2-3 decades 
in Ml I Mi, and AMjMi. 

We also investigated the dependence of the conditional 
mass functions on the choice of group finding algorithm. 
Using the "best fit" h^a values from §5, the analytical distri- 



butions fit about as well for FOF(0.3) as for FOF(0.2), while 
for FOF(0.15) and SO(180), the fits were slightly worse. We 
show the results for the n=-l spherical overdensity groups 
in Figure 7, with h^a = 1.96. In this case, the analytical 
distribution actually fits even better if 8co = 1.69 is chosen. 



7 FORMATION TIMES 

In hierarchical models halos evolve continuously by accret- 
ing smaller halos and by merging with comparable and larger 
halos. Thus there is no clear cut way of defining when a par- 
ticular halo formed. As a working definition we have adopted 
the formation time to be the point at which half the mass 
of a halo is assembled. The distribution dpi/dti{ti\M2,t2), 
equation (2.19), gives the probability that half the mass of 
a halo of mass M2 identified at time t2 was assembled at an 
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Figure 6. Like Figure 4, but for the n = model. The values of Aa/ai indicated on the figures again correspond to the intervals over 
which increases by factors \/2 and 16 respectively. 



earlier time tf. Here we compare this analytical distribution 
of formation times with estimates made directly from our 
N-body simulations. 

Starting with a set of group catalogues defined using one 
of the group finding schemes of §5 we proceed as follows. For 
each group of mass M identified at a particular output epoch 
at which the expansion factor equals ao , we identify its most 
massive progenitor at all earlier output times. Progenitors 
of a group are defined to be all those groups present at an 
earlier epoch that have the majority (normally greater than 
90%) of their particles incorporated into the final group. We 
then determine at which epoch the mass of this most mas- 
sive progenitor first becomes larger than M/2. This epoch 
is taken to define the formation time of the halo and de- 
fines the expansion factor at at formation. The histogram 
of af values built up for a particular choice of ao and M 
defines the formation time distribution dpi /din ai\M,ao - For 



our scale free simulations, this distribution is a function of 
only M / and af/ao. Hence we can once again combine 
the histograms from different final output times, ao, to bet- 
ter determine the numerical estimate of dpi /din ai\M/Mi^- 
We choose to use only final epochs separated by a factor 2 
in M^, so that they are approximately independent, and to 
estimate Poisson errors from the combined number of halos 
in each bin. 

Figures 8,9 and 10 compare these results for FOF(0.2) 
groups for a range of M/M^ with the analytical predictions 
of equation (2.19), for spectral indices n = —2,-1 and 
respectively. In constructing these plots, we use only halos 
with TV > 40 particles at the epoch ao. (Note, if the epoch 
ao is identified with the present epoch then the quantity 
plotted along the x-axis, log(af/ao) = — log(l -|- Zf), where 
Zi is the formation redshift). In each case, there is a clear 
trend for larger M / M^^ halos to be both younger on average 
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Figure 7. Like Figure 5, but for SO(180) groups and with the best fit value of Sco = 1.96 adopted from Figure 2. 



and to have a narrower range of ages. This behaviour, and in 
fact the precise shape of the distributions, is well reproduced 
by the analytical distributions given by equation (2.19). We 
find agreement over 2-3 decades in M/M^. 

We also investigated the formation times of groups iden- 
tified using the FOF algorithm with differing values of the 
linking length b and for groups defined by SO with k = 180. 
Once the threshold 8co was adjusted to the appropriate "best 
fit value" found for the mass functions, then the agreement 
between the analytical and numerical distributions was as 
good as for the FOF(0.2) groups. In contrast to the fits to the 
conditional mass functions in §6, the fits to the formation 
time distributions were appreciably worse for FOF(0.15), 
FOF(0.3) and SO(180) groups if the value Sco = 1.69 was 
used instead of the appropriate "best fit" value. As an ex- 
ample we show the distributions of formation times for the 
case of the SO(180) groups in Figure 11. 



In paper I, we defined a scaled variable 
_ ge(a/) - ^e(ao) 

(M/M4ao))("+^)/^ (ao/gf-l) 

(2(" + 3)/3 _ 1)1/2 

(the second line is for self-similar models) in terms of 
which the analytical distribution of formation times dpi I dCoi 
is independent of M/M^ and very nearly independent of 
the spectral index n. This last property was demonstrated 
graphically in Figure 7 of Paper I. Here we transform each 
of distributions from Figures 8,9 and 10 into distributions 
dpi I dCoi which we show in Fig. 12. These figures confirm 
that for the groups found in the N-body simulations there 
is also very little dependence of these distributions on either 
mass or spectral index. The short-dashed curves in Fig. 12 
are the analytical distributions while the long-dashed curves 
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Figure 8. Comparison of the distribution of formation epochs derived from the FOF(0.2) groups in the n = —2 simulation (sohd curves) 
with the analytical prediction (2.19) (dashed curves), for the various values of M/M^ indicated. We have assumed Sco = 1.69. For the 
N-body curves, halos were identified and formation epochs a ^ found for a set of identification epochs ao differing by powers of 2 in Af^, 
and the results averaged. Only halos with > 40 were used. The error bars represent Poisson errors corresponding to the total number 
of halos in each bin in aj^, after combining the different epochs. 




Figure 9. Same as Figure 8, but for the n = 



— 1 model. 
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Figure 10. Same as Figure 8, but for the n = model. 
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Figure 11. Same as Figure 9, but for SO(180) groups in the n = 



— 1 model, and with Sco = 1-96. 
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Figure 12. Distributions of dpi / dCoi for the N-body simulations compared with the analytical and Monte Carlo predictions from Paper L 
The thin solid lines show the N-body results for different values of M/M^ (averaged over output times as in Figures 8 - 10). The N-body 
curves plotted are for the mass ranges M/M^ = 0.9 — 120, 0.1 — 15, 0.2 — 7 for n = —2, —1, respectively. The short dash and long dash 
curves are respectively the analytical and Monte Carlo predictions. 



are the Monte-Carlo distributions taken from Paper L The 
Monte Carlo method attempts to follow a single path back 
through the merger tree leading to a given object, choos- 
ing the more massive progenitor each time a halo splits 
in two, back to when the largest progenitor has fallen be- 
low half the final mass. Clearly, the distributions estimated 
from the N-body simulations are well described by the an- 
alytical predictions (equation 2.19), while the Monte Carlo 
model tends to overestimate halo ages. The problem with 
the Monte Carlo method appears to lie with an incorrect 
weighting of the probabilities for the distribution of pro- 
genitor masses at each branching point. We are currently 
working on an improved Monte Carlo procedure which gives 
results very close to the analytical and N-body ones, and 
this will be described in a future paper. 



8 CONCLUSIONS 

We have tested the statistical predictions of the Press- 
Schechter model for the mass function of dark matter halos 
(Press & Schechter 1974; Bond et al. 1991; Bower 1991), and 
its extension to halo lifetimes and merger rates (Paper I), 
against a set of large N-body simulations. The simulations 
used 128'^ particles and modelled self-similar clustering with 
Q = 1 and initial power spectra P{k) oc , with spec- 
tral indices n = —2,-1,0. The comparison reveals that in 
every respect the analytical formulae produce remarkably 
good fits to the numerical results. Although tested for self- 
similar clustering, the analytical formulae are also expected 
to apply for arbitrary Q and more general power spectra, 
provided that structure grows hierarchically from Gaussian 
density fluctuations in cold coUisionless matter. 

Dark matter halos were identified in the simulations 
using two alternative methods. The first was the standard 



percolation or "friends of friends" method, which effectively 
selects objects bounded by surfaces of specified density. We 
investigated linking lengths b = 0.15, 0.2 and 0.3 in units 
of the mean interparticle separation, corresponding to mean 
halo overdensities in the range ~ 50 — 500, smaller b corre- 
sponding to higher density. The second, "spherical overden- 
sity", method finds spheres of a specified mean density k, 
where we used k = 180. 

To apply the analytical formulae for mass functions and 
merging, three choices have to be made: First one must se- 
lect the form of spatial filter which is applied to the linear 
density field in order to define cr{R), the r.m.s. fluctuation 
as a function of length scale. Secondly, one must choose a 
relation between filter radius R and mass M to derive a{M) 
from a{R). Finally, one must set the critical density thresh- 
old for collapse, 8co- We have investigated top hat, sharp 
A;-space and Gaussian filtering, with a mass-radius relation 
M = nJjR'' for some filter-dependent constant, 7f. For top 
hat and sharp A;-space filtering, the value of ^co required 
to best match the analytical Press-Schechter mass function 
with the N-body results is independent of the spectral index 
of the linear density field for —2 < ra < 0, when n is taken 
to be the value obtained by integrating the window function 
over all space (for the top hat, M is just the mass enclosed 
within radius R). For Gaussian filtering, on the other hand, 
this independence holds only if n is taken to be a factor 
~ 2.5 larger than is given by this integral, contrary to what 
has been assumed in previous work. When halos are selected 
using percolation with b = 0.2, which selects halos having 
mean overdensity ~ 100 — 200, the best fitting 8co values 
for each of the filters are within 20% of the value 8co = 1.69 
predicted by the analytical model for the collapse of a spher- 
ically symmetric overdense region in an Q = 1 universe. In 
this case, the Press-Schechter mass function fits the N-body 
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results to an accuracy of ~ 30%. When the group selec- 
tion method is changed, the best fit 8co also changes. The 
mass functions for the b = 0.15 and k = 180 groups require 
larger values of 8co (further from 8co = 1.69) in the Press- 
Schechter formula to give reasonable fits, compared to the 
b = 0.2 groups, and even then fit somewhat less well. 

Overall, there seems little to choose between the dif- 
ferent types of filtering, but we have concentrated on the 
top hat filter in most of our comparisons because this is 
more standard. For the conventional choices of top hat fil- 
tering with 8co = 1.69, the analytical mass functions differ 
by less than a factor 1.5 — 2 from those estimated from the 
b = 0.2 percolation groups in the simulations, over a range 
of 10'^ in mass (see Fig. 1). The error is largest for the rare 
high-mass halos, but is typically < 30% for the more nu- 
merous halos which contain most of the mass. The error in 
the mass function can in fact be reduced to < 30% overall 
for this case by increasing 8co by ~ 10%. The comparison 
is limited to the abovementioned range of mass by the dy- 
namic range of our simulations, which span a factor 10'^ in 
mass from the smallest resolved halos (containing at least 
20 particles) to the most massive. Although still limited in 
dynamic range, this is a considerable improvement over the 
comparison made by Efstathiou et al. (1988), which utilized 
simulations containing only l/64th of the number of parti- 
cles of our simulations. 

With the same choices of filter, threshold, and group- 
finding, we also find remarkable agreement between the halo 
merger rates measured from the simulations and the analyt- 
ical predictions. All the trends seen in the dependence of the 
merger rate on the masses of the two halos involved and on 
the time interval considered are reproduced quantitatively 
by the analytic formula (2.17) (see Figs. 4-6). In fact, equa- 
tion (2.17) is a reasonable fit to the numerical estimates 
over the full range of masses, roughly 2-3 decades, that we 
are able to explore. A similarly impressive agreement, again 
for a wide range of masses, is seen when one compares the 
distribution of halo formation times estimated from the sim- 
ulations with the analytical formula (2.19) (see Figs. 8-10). 
For the b = 0.15, b = 0.3 and k = 180 groups, the agreement 
for merger rates and formation times is about as good as for 
the b = 0.2 groups, provided that instead of the standard 
value 8co = 1.69, one uses the values of 8co which best fit 
the mass functions in the other formulae too. 

We end with a caveat. Despite its success in matching 
the results of N-body simulations, the Press-Schechter ap- 
proach from which our formulae are derived falls some way 
short of being a rigorous analytical model of gravitational 
instability and non-linear dynamics. It is instead based on 
the ansatz (Bond et al. 1991) that the mass in non-linear 
objects of mass M can be equated with the mass within re- 
gions whose linear theory density perturbation 8 exceeds a 
threshold value, 8c, when smoothed on the mass scale M, 
but is below the threshold on all larger scales. In particu- 
lar, no regard is paid to the shape or size of these regions. 
Many will enclose less than mass M, making it impossi- 
ble that they form objects of mass M without combining 
with other nearby material. It is clear that any physically 
realistic model must take account of the properties of the 
linear density field across the entirety of each region that 
collapses to form a non-linear halo. The "peak-patch" anal- 
ysis of (Bond & Myers 1993a) is the first model that ad- 



dresses this aspect of the problem. The disadvantage of this 
more rigorous treatment is that it is very complex and does 
not lead to analytical formulae for halo mass functions or 
merger rates. Thus, despite the fundamental flaws in the 
Press-Schechter approach, the analytical formulae that it 
yields are extremely valuable if they are an accurate descrip- 
tion of the true halo mass functions and merger statistics. 
This work confirms that these statistical predictions repro- 
duce remarkably well the non-linear hierarchical evolution 
of dark matter halos in large scale-free cosmological N-body 
simulations. They therefore provide a valid and extremely 
useful framework in which to study galaxy and cluster for- 
mation in hierarchical models. 
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